%function runCase ()
%
%
clear all; clc;

datafile = 'case162.m'; 
sim = 1;               % 0: reuse older data on initial shedding time
endTime = 1.8;

% Run initial.m & Ymatrix.m
initialize( datafile );
Ymatrix();             % Y3fun: function to change Y by cutting rate

% Simulate before cutting time
simulate(sim, endTime);

tic();
% Model and initialize Optimal Control Problem
modelOCP( endTime );
preprocessOCP();       % NLP options, guess options

% Optimazation
solveNLP(); 
toc();
%% Get results
global dataSim postSimFunc        % simulation results
global sol                        % solution of "case$x.m"
% Set ideal uG and uL; Set No. of generator to display
% sol.uG = []; sol.uL = [];
idx = (1 : length(sol.uG))'; nidx = length(idx);            % --- choose "delta" to display
% Whole time without cutting gen & load
deltaRec = []; load(dataSim, 'tRec', 'deltaRec', 'coiRec');
cv = deltaRec(:,idx) - repmat(coiRec, 1, nidx); % curve
% Solution of NLP: time phase 3
t3 = linspace(sol.time(1), sol.time(end), 50)';
[x3, coi3] = interpAt(t3, sol, 'states', idx);
cv3 = x3 - repmat(coi3, 1, nidx);
% Post simulation varification
[tP, dltP, coiP] = feval(postSimFunc, sol.uG, sol.uL);
cvP = dltP(:,idx) - repmat(coiP, 1, nidx);
% Plot and Print
plot(tRec, cv, t3, cv3, tP, cvP);
axis([0, t3(end), -150, 150]);
fprintf('Gen Cut:\n'); fprintf('%8.4f', sol.uG'); fprintf('\n');
fprintf('Load Cut:\n'); fprintf('%8.4f', sol.uL'); fprintf('\n');

%end